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Abstract: The freshwater gudgeon Huigobio chenhsienensis (Cypriniformes: Gobionidae) is a small fish 
endemic to southern China. In this study, we used mitochondrial cytochrome b gene (Cytb), from 
wide-ranging samplings of H. chenhsienensis from the Ou River (the central of southern China) to the 
Yangtze River Basin (the northernmost part of southern China) to explore genetic variations and the 
evolutionary history of H. chenhsienensis in southern China. In total, 66 haplotypes were identified from 
Cytb sequences of 142 H. chenhsienensis individuals, which could be divided into lineages A, B, and C 
with divergence times of ~4.24 Ma and ~3.03 Ma. Lineage A was distributed in the lower reaches of the 
Yangtze River, the Oujiang River, and the Jiao River, lineage B was distributed in the Qiantang River and 
the Cao’e River, whereas lineage C was restricted to the Poyang Lake drainage from the middle reaches 
of the Yangtze River. Lineage A could be subdivided into sub-lineages A-I, A-II, A-III, and A-IV, with 
divergence times of 1.30, 0.97, and 0.44 Ma. Lineage C could be subdivided into sub-lineages C-I and 
C-IIL, with a divergence time of 0.85 Ma. Our findings indicate that climate change during the Pliocene 
and Pleistocene eras, as well as the limited dispersal ability of H. chenhsienensis, have been major drivers 
for shaping the phylogeographical patterns of H. chenhsienensis. 
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1. Introduction 


Phylogeography focuses on the geographic distributions of inter- and intraspecific 
lineages, as well as the historical processes that have affected their evolution [1,2]. Fresh- 
water fish in contemporary isolated drainages are often the products of historical river 
evolutions. Therefore, phylogeographical studies of freshwater fish are significant for 
understanding the historical reorganization or rearrangements of rivers within and among 
drainages [3-6]. The genetic diversity observed within a species is a key facet of biodiver- 
sity that reflects the evolutionary history and adaptability of a species as a whole, or its 
constituent populations [7,8]. 

Grant and Bowen [9] classified four categories of genetic diversity patterns according 
to levels of polymorphism for haplotypes and nucleotides. They defined high or low 
haplotype polymorphism (h) as h > 0.5 or h < 0.5, and high or low nucleotide polymorphism 
(zt) as 7¢ > 0.005 or 7t < 0.005. The Yangtze River is the northern boundary for freshwater fish 
distributed in southern China [10]. Previous studies (Table S1) have shown two categories 
of genetic diversity for freshwater fish in southern China, i.e., type II (high value and low 
7t value) and type IV (high h value and high z value), as defined by Grant and Bowen [9]. 
The potential formation mechanism of type II cases is a genetic bottleneck followed by rapid 
population expansion [9]. However, widespread climate oscillations in the Pleistocene and 
secondary contacts resulting from population expansions after the glacial period may have 
shaped the type IV genetic diversity pattern [11]. 
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Based on the geographical distribution of lineages and the degree of genetic differenti- 
ation, Avise [12] summarized five categories of phylogeographic patterns as follows: “Deep 
gene tree, major lineages allopatric” (pattern I), “Deep gene tree, major lineages broadly 
sympatric” (pattern II), “Shallow gene tree, lineages allopatric” (pattern III), “Shallow gene 
tree, lineages sympatric” (pattern IV), and “Shallow gene tree, lineage distributions varied” 
(pattern V). Prior studies on freshwater fish in southern China (Table S2) have shown 
three categories of phylogeographic patterns. The main driver of pattern I is the long-term 
isolation of drainages resulting from biogeographic barriers formed by geological tectonic 
movement [13-15]. The main driver of pattern II is river capture or historical connections 
of the river facilitated by sea level fall during the Pleistocene glacial period [16,17]. Pattern 
V is relatively rare in nature, and is attributed to the low contemporary gene flow between 
populations which have historically been closely connected [12,18]. Although Chinese 
ichthyologists have carried out phylogeographical studies on a variety of fish in southern 
drainages, it is still necessary to add phylogeographical cases of fish species in the above 
area so that more evidence can be accumulated to understand the evolutionary processes 
that generate and maintain intraspecific genetic variation, and lay a scientific foundation 
for the effective management and rational protection of genetic resources for freshwater 
fish in China. 

Huigobio chenhsienensis Fang 1938 is classified under the order Cypriniformes, family 
Gobionidae, subfamily Gobioninae, and genus Huigobio [19,20]. It is a species endemic 
to southern China, and is distributed from the Ou River to the Yangtze River [20-23]. 
H. chenhsienensis is small (standard length of less than 8 cm) and usually inhabits the 
upper reaches of streams, feeding mainly on benthic algae [21,22,24]. The small body size 
and stream headwater habitats may indicate the potential limited dispersal ability of H. 
chenhsienensis. In addition to H. chenhsienensis, the genus Huigobio comprises two other 
species: H. exilicauda Jiang and Zhang 2013 and H. heterocheilus Sun, Li, Tang and Zhao 
2022. H. exilicauda is distributed in the Pearl River, whereas H. heterocheilus is found in the 
Xiang River, a tributary of Dongting Lake drainage in the middle reaches of the Yangtze 
River [22,23]. 

In this study, mitochondrial cytochrome b gene (Cytb) sequences were obtained from 
samples covering the distribution drainages of H. chenhsienensis. Based on analyses of 
phylogeny, genetic structure, and population history, we discuss the roles of river evolution 
and paleoclimatic events in the formation of phylogeographic patterns of H. chenhsienensis 
in southern China. 


2. Materials and Methods 
2.1. Sample Collection 


From November 2011 to September 2019, 142 specimens of H. chenhsienensis were 
collected from the Ou River, Jiao River, Qiantang River, Cao’e River, and Yangtze River, 
using small set nets and gill nets by ourselves or with help from local fishermen. There 
are a total of 22 localities, with 1 to 18 specimens for each locality (Figure 1 and Table 1). 
In order to take into account fish euthanasia, the experimental fish were anesthetized to 
lose consciousness using the Eugenol. First, the Eugenol (purity > 99.99%) was diluted 
into an aqueous solution with a concentration of 0.25 mL-L7!. Second, experimental fish 
were placed in the aqueous solution of the Eugenol until they were anesthetized to lose 
consciousness. Third, the anesthetized fish were fixed in 75% and then transferred to 95% 
ethanol for long-term storage. All specimens collected in this study were deposited in the 
Zoological Museum of Fudan University. 


2.2. Sequence Obtaining and Processing 

Genomic DNA was extracted from muscle tissue using a salting-out extraction 
method [25]. A pair of primers (CGluF: 5’/-TCTTGCTCGGAYTTTAACCGAG-3’, 
CThrR: 5‘-TCTTGCTCGGAYTTTAACCGAG-3’) [26] were used for PCR amplification 
of Cytb gene. The PCR amplification was conducted in a 30 uL reaction volume which 
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contained 1.2 uL of template DNA, 0.6 uL of each primers, 15.0 wL of 2 x Es Taq MasterMix, 
and 12.6 uL of ddH20. The reaction conditions of PCR amplification for Cytb gene were the 
same as Chai and Fu [27] except that the primers were annealed at 58 °C. Moreover, the 
primers and annealing temperature for PCR amplification of the mitochondrial genome 
used in this study are listed in Table S53. The PCR products were subsequently purified and 
sequenced by Jieli Biology Co., Ltd., Shanghai, China. 
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Figure 1. Localities of Huigobio chenhsienensis specimens and spatial distributions of mitochondrial 
cytochrome b gene (Cytb) lineages. The green (paleo-Yangtze River), magenta (paleo-Jiao River), and 
orange (paleo-Ou River) lines represent paleo-channels of coastal rivers’ downward movement to 
continental shelf at an elevation of —130-0 m during the last glacial maximum (LGM) when the sea 
level fell by about 130 m. The code of each population is shown in Table 1 and lineages are defined in 
Figure 2. FJ—Fujian Province; JX—Jiangxi Province; ZJ—Zhejiang Province; AH—Anhui Province. 


The obtained sequences were assembled in Sequencher v5.4 (Gene Codes, Ann Arbor, 
MI, USA) and supplemented with manual correction for incorrect base recognition of the 
chromatograms. Sequence alignment was performed with the default algorithm under 
MAFFT v7.427 [28]. Amino acid translation was performed in DAMBE v7.2.43 to check 
whether there were terminators in the sequences [29]. The newly generated Cytb sequences 
were imported into DnaSP v6.12.01 to obtain haplotype data [30]. Basic information on 
the Cytb sequence was acquired using MEGA (Molecular Evolutionary Genetics Analysis) 
version 7.0.26 (Arizona State University, Philadelphia, PA, USA) [31]. 
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Figure 2. A time-calibrated phylogeny and ancestral area reconstructions among 66 Cytb haplotypes 
of H. chenhsienensis. The mean value of the divergence times and their 95% confidence interval are 
shown above the nodes. The posterior probabilities are below the nodes. The numbers labelled in 
black circles represent the serial numbers of the major nodes (posterior probability more than 0.9). 
The outgroup taxon (Huigobio heterocheilus) is not shown, and its information in detail is provided in 
Table S4. The color block at each node reflects ancestral distribution. 
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Table 1. Information on localities of specimens and mitochondrial cytochrome b gene (Cytb) haplo- 
types for Huigobio chenhsiensis. 


Basin and Locality Lat. Long. N Cytb Haplotype 

Ou River 

1. Yongjia County, Zhejiang Province 28.35° N 120.75° E 14 BN BO): rand He a BOGE) BOM), 
2. Jinyun County, Zhejiang Province 28.66° N 120.09° E 2 A11(2) 

Jiao River 

3. Tiantai County, Zhejiang Province 29.14° N 121.03° E 1 A12(1) 

Cao’e River 

4. Xinchang County, Zhejiang Province 29.51° N 120.91° E 15 BO1(4), BO2(1), BO3(8), BO4(2) 

Qiantang River 

5. Pan’an County, Zhejiang Province 29.06° N 120.45° E 1 BO5(1) 

6. Yugian Town, Lin’an City, Zhejiang Province —30.19° N 119.40° E 6 B06(2), BO7(1), BO8(2), BO9(1) 

> oe Hover Taam aly ehepans 30.177N 119.22°E —-10 B08(3), B10(4), B11(1), B12(1), B13(1) 

8. Yunti Village, Ningguo City, Anhui Province —30.37° N 119.21° E 11 B05(1), BO8(2), B12(1), B14(1), B15(1), B16(4), B17(1) 
9. Huangshan City, Anhui Province 29.72° N 118.31° E 1 BO5(1) 

10. Jiangshan City, Zhejiang Province 28.76° N 118.66° E 1 B18(1) 

Yangtze River 

11. Anji County, Zhejiang Province 30.64° N 119.67° E 6 A13(5), A14(1) 

12. Ningguo City, Anhui Province 30.62° N 118.95° E 7 A15(1), A16(1), A17(1), A18(1), A19(1), A20(1), A21(1) 
13. Shitai County, Anhui Province 30.21° N 117.50° E 10 cae ace a eee oan sil Ge ac eal 
14. Fuliang County, Jiangxi Province 29.52° N 117.28° E 1 C01(1) 

15. Yushan County, Jiangxi Province 28.69° N 118.28° E 1 C02(1) 

16. Fengxin County, Jiangxi Province 28.71° N 115.38° E 6 C03(1), C04(2), C05(1), C06(1), CO7(1) 

17. Jing’an County, Jiangxi Province 28.85° N 115.39° E 18 C04(2), C08(13), CO9(1), C10(1), C11(1) 

18. Tonggu County, Jiangxi Province 28.53° N 114.39° E 9 C06(5), C12(2), C13(1), C14(1) 

19. Xiushui County, Jiangxi Province 29.02° N 114.54° E 1 C15(1) 

20. Ningdu County, Jiangxi Province 26.46° N 116.00° E 5 C16(3), C17(2) 

21. Suichuan County, Jiangxi Province 26.31° N 114.52° E 15 C16(9), C17(2), C18(1), C19(3) 

22. Anfu County, Jiangxi Province 27.39° N 114.63° E 1 C16(1) 


Lat.—latitude; Long.—longitude; N—number of individuals; number of individuals for each haplotype is 
in bracket. 


2.3. Phylogeographic Analysis and Divergence Time Estimation 


The Cytb time tree of H. chenhsienensis was reconstructed in BEAST v2.6.6 software 
(University of Auckland, Auckland, New Zealand) [32] based on haplotype sequences. 
Taking full account of the lack of fossil information on Huigobio for estimating the differen- 
tiation time of H. chenhsienensis, a mitochondrial genome time tree of Huigobio and its close 
relatives was thus developed to obtain the differentiation time of H. chenhsienensis and H. 
heterocheilus as the secondary calibration point by using two fossil species of Gobionidae 
with complete skeletons and corresponding stratigraphic age information (data sources 
and results are shown in Table S4 and Figure $1). The best base substitution model for the 
13 mitochondrial protein genes was obtained based on the automatic search function of 
the bModelTest module in BEAST v2.6.6, the Birth-Death model was selected, and a tree 
model was set at the same time [33]. The first fossil species corresponded to the divergence 
node of Gobionidae and Acheilognathidae [34]. This fossil species, which occurred in the 
Middle Eocene (50.5-42.0 Ma) [35] is tPalaeogobio zhongyuanensi [36], i.e., calibration point 1 
(C1), and was set to LogNormal distribution, 1. = 3.82, o = 0.05. The second fossil species 
corresponded to the divergence node of Gnathopogon and Coreoleuciscus [37]. This fossil 
species occurred in the Miocene (20.4-15.0 Ma) [38] and is tGnathopogon macrocephala [36], 
i.e., calibration point 2 (C2), and was set to LogNormal distribution, 1 = 2.86, o = 0.078. 
The molecular clock was set as a relaxed clock. Two independent runs of 400 million 
generations with sampling every 10,000 generations in each run were performed in the 
Bayesian analysis. The two runs were combined and the first 30% sampling trees were 
discarded as burn-in. Effective sample size (ESS) was evaluated using Tracer v1.7.0 [39]. 
When the ESS was greater than 200 for each parameter, the results were considered to 
be convergent. 
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To test for the geographic correspondence of haplotypes across all populations of H. 
chenhsienensis, the median-joining network was constructed under NETWORK v10.2.0.0 
(Fluxus Technology Ltd., Colchester, England; https: / /www.fluxus-engineering.com, ac- 
cessed on 15 January 2022). Combined with the geographical distribution and individual 
number of haplotypes, it was further modified in Adobe Illustrator. 

The interspecific divergence time between H. chenhsienensis and H. heterocheilus was 
11.66 Ma (95% HDP: 13.04—10.44 Ma), which was obtained from the time tree of the mito- 
chondrial genome presented above and used as a secondary calibration point to infer the 
Cytb time tree of H. chenhsienensis. The specific settings were as follows: normal distribu- 
tion, m = 11.74, s = 0.66. Using the likelihood-ratio test in DAMBE v7.2.43, the result of 
the molecular clock test (X* = 96.11, df. = 69, p = 0.0172) did not support the molecular 
clock model [29], and the relaxed-clock model was set. H. heterocheilus was selected as 
an outgroup. The Birth-Death model was selected as a tree prior. Two independent runs 
of 100 million generations, with sampling every 1000 generations, were performed. The 
remaining settings were the same as the time tree construction above. 


2.4. Genetic Structure and Population History 


The genetic diversity of H. chenhsienensis was estimated, and the genetic differentiation co- 
efficients, Nst and Ggz, of the populations were calculated and compared with determination 
of the presence of a significant genetic structure using DnaSP v6.12.01 [30,40]. The total Bsr 
and pairwise comparison Psy values among basins were obtained in Arlequin v3.5.1 [41]. To 
detect the spatial genetic structure of H. chenhsienensis populations, the maximum ®cy7 value 
and its corresponding grouping were obtained by simulating the groupings of different popu- 
lations using SAMOVA v2.0 [42]. For each lineage, Arlequin v3.5.1 was applied to conduct 
neutrality tests [43,44] and mismatch distribution tests with a spatial expansion model [45]. 
Bayesian skyline plot analysis of each lineage was implemented in BEAST v2.6.6 [32], the 
Relaxed Clock Log-Normal distribution model was selected, and the tree model was set as 
the Coalescent Bayesian Skyline. The best base substitution model used for the analysis was 
obtained from jModeltest v2 [46]. The Cytb mutation rate (1.30% substitutions per site per 
million years) obtained from the time tree of Cytb haplotypes of H. chenhsienensis was used in 
the reconstruction of the Bayesian skyline plot. 


2.5. Ancestral Area Reconstruction 


The ancestral area reconstruction of H. chenhsienensis was performed by loading the 
BioGeoBEARS package version 0.2.1 [47] in the R v3.5.0 program. The analysis used the 
Cytb haplotypes time-calibrated tree of H. chenhsienensis as the input tree and the basins 
as the geographic distribution unit. Finally, the optimal model was selected according to 
Akaike information criterion (AIC) to reconstruct the ancestral area of H. chenhsienensis [48]. 


2.6. Inferring Palaeodrainage along the Coast 


Topographic and bathymetric information on a 30 arc-second interval grid at a global 
scale, with accuracy in meters, was retrieved from the digital elevation model (DEM) of the 
General Bathymetric Charge of the Oceans 2021 Grid (GEBCO_2021 Grid; http://www. 
gebco.net/, accessed on 10 November 2021). The maximum sea level decline during the 
glacial period of the Late Pleistocene was about —130 m [49,50]. Therefore, the elevation 
and bathymetry data of the maximum sea level retreat during the glacial period of the 
Late Pleistocene could be used to infer potential connections between palaeodrainages and 
current basins along the coast of China through the Arc Hydro Tools Box in ArcGIS v.10. The 
palaeodrainage reconstruction methods were performed in reference to Thomaz et al. [51] 
with the following steps: first, the DEM of GEBCO_2021 Grid was imported into ArcGIS 
v.10; secondly, the —130 m basic contours were drawn using the Contour List tool in the 
Surface Tools box to estimate the maximum extent of continental shelf exposure during the 
Pleistocene, followed by the generation of a clipped digital elevation model grid chart using 
the Fill, Flow Direction, Watershed, Flow Accumulation, and Stream Order tools in the Arc 
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Hydro Tools Box. Finally, the raster streams were converted to vectors, and a visualized 
palaeodrainage reconstruction map was generated using the Stream to Feature tool. 


3. Results 
3.1. Phylogenetic Relationships and Divergence Time 


A total of 66 haplotypes (Table 1, GenBank Accession: OM647979-OM648044) were de- 
tected among the Cytb sequences of 142 H. chenhsienensis specimens. All the Cytb sequences 
had an identical length (1140 bp), with a GC content of 41.3%, containing 184 variable sites 
and 161 parsimony informative sites. The time-calibrated phylogeny (Figure 2) showed that 
66 Cytb haplotypes of H. chenhsienensis could be divided into lineages A (29 haplotypes), 
B (18 haplotypes), and C (19 haplotypes); divergence times among them were 4.24 Ma 
and 3.03 Ma. The lineage A could be subdivided into four sub-lineages, namely, A-I, A-IL, 
A-III, and A-IV, which diverged from each other at 1.30, 0.97, and 0.44 Ma. The lineage C 
consisted of two sub-lineages, C-I and C-II, with a divergence time of 0.85 Ma. 

There were 61 mutational steps between lineages A and B + C and 54 mutational steps 
between lineages B and C observed in the haplotype network. All 66 haplotypes were 
private haplotypes for each river basin (Figure 3). In lineage A, the three sub-lineages, A-I, 
A-II, and A-III, were only distributed in the lower reaches of the Yangtze River, whereas 
the sub-lineage A-IV was distributed across the lower reaches of the Yangtze River, Ou 
River, and Jiao River (Figure 3a). The lineage B was distributed in the Qiantang River 
and Cao’e River (Figure 3b). The lineage C was only distributed in Poyang Lake of the 
middle reaches of the Yangtze River (Figures 1 and 3c), among which the sub-lineage C-I 
was distributed in Gan River of the Poyang Lake river system (Locality No. 20-22, as 
shown in Table 1 and Figure 1) and the sub-lineage C-II was distributed in Rao River, Xin 
River, and Xiu River of the Poyang Lake river system (Locality No. 14-19, as shown in 
Table 1 and Figure 1). 


3.2. Genetic Structure 


Among the distribution basins of H. chenhsienensis (excluding Jiao River with one 
specimen), haplotype diversity (Table 2) was highest in the Ou River (i = 0.9500) and lowest 
in Cao’e River (i = 0.6667); nucleotide diversity (Table 2) was highest in the Yangtze River 
(zt = 0.0371) and lowest in Cao’e River (7t = 0.0007). 

Comparison between two genetic differentiation coefficients (Nst > Gsr: Nest = 0.796, 
Gsr = 0.072) indicated a significant geographic structure for H. chenhsienensis. The total ®sr 
among basins was 0.627 (p = 0.000), and the pairwise comparison Psy values among basins 
were relatively high and significant (Table 3). The SAMOVA divided the five basins into 
four genetic groups as follows: Yangtze River, Qiantang River, Cao’e River, Ou River + Jiao 
River, where ®cz value (=0.741) was the largest. 


3.3. Demographic History 


Tajima’s D value and Fu’s Fs value of lineage A for H. chenhsienensis were insignifi- 
cantly positive and marginally significantly negative, respectively, and both indices were 
significantly negative in lineage B and insignificantly positive in lineage C (Table 4). Bi- 
modal distributions were observed in lineages A and C, whereas lineage B exhibited a 
unimodal peak (Figure 4a). The Bayesian skyline plot analysis (Figure 4b) revealed that the 
population expansion of lineages A and B started at 0.088 Ma and 0.110 Ma, respectively; 
however, lineage C contracted from 0.040 Ma. 


3.4. Ancestral Area Estimation 


By evaluating six biogeographic models of BioGeoBEARS analyses, the DIVALIKE+J 
model provided an optimal fit to the data based on AIC scores (Table S5). The ancestral area 
reconstructions based on the optimal model supported the suggestion that the common 
ancestor of H. chenhsienensis was widely distributed in the middle reaches of the Yangtze 
River, the lower reaches of the Yangtze River, the Qiantang River, and the Cao’e River at 
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4.24 Ma (95% HPD: 6.05-2.62 Ma). The current geographic distribution pattern was formed 
after three vicariant and three dispersal events (Figure 2). 
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Figure 3. A haplotype network of lineages A (a), B (b) and C (c) for H. chenhsienensis. Missing 
haplotypes are represented with black dots. The number of mutational steps is marked between the 
haplotypes unless that number is one. 


Table 2. Genetic diversity of H. chenhsienensis. 


Basin Nisnbes Number of Number of Private Haplotype Nucleotide 
HaploTypes (Nh) Haplotypes (Nph) Diversity (h) Diverstiy (7z) 
Ou River 16 11 11 0.9500 + 0.0364 0.0028 + 0.0017 
Jiao River 1 1 1 - - 

Cao’e River 15 4 4 0.6667 + 0.0991 0.0007 + 0.0006 
Qiantang River 30 14 14 0.9126 + 0.0307 0.0024 + 0.0015 
Yangtze River 80 36 36 0.9370 + 0.0149 0.0371 + 0.0180 
Total 142 66 66 0.9722 + 0.0054 0.0460 + 0.0222 
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Table 3. Pairwise comparison Psy values among basins (lower diagonal) and associated p values 
(upper diagonal) of H. chenhsienensis. 


Ou River Cao’e River Qiantang River Yangtze River 
Ou River 0.000 0.000 0.000 
Cao’e River 0.986 0.000 0.000 
Qiantang River 0.981 0.367 0.000 
Yangtze River 0.546 0.510 0.559 


Jiao River was excluded from the analysis because it had only one sample. Significant p values (p < 0.05) by 
Bonferroni correction are displayed in bold. 


Table 4. Neutrality tests for lineages A, B, and C of H. chenhsienensis. 


Lineage Tajima’s D Fu’s Fs 
A 0.320 (p = 0.706) —4.231 (p = 0.095) 
B —1.768 (p = 0.014) —9.126 (p = 0.000) 
Cc 0.528 (p = 0.781) 1.220 (p = 0.714) 
(b) 
Lineage A SSD: 0.023 (p= 0.198) - Lineage A 


Rag. index: 0.009 (p =0.987) 
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Figure 4. Population history including mismatch distributions (a) and Bayesian skyline plots (b) for 
lineages A, B, and C of H. chenhsienensis. In the mismatch distributions, the bar graphs represent 
observed values, and the black solid dots and solid lines represent expected values under the 
population expansion model. For the Bayesian skyline plots, the actual value of effective population 
size is represented with the solid blue line, 95% confidence interval is represented with the shaded 
range, and the time range of population expansion or recession is shown with the red box. 
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3.5. Palaeodrainage Reconstructions 


Among the distribution basin of H. chenhsienensis, the reconstruction of the coastal 
palaeodrainages (Figure 1) showed that the Yangtze River, Qiantang River, and Cao’e River 
were interconnected on the exposed continental shelf during the Last Glacial Maximum, 
whereas the Ou River and Jiao River entered the sea independently during the same period. 


4, Discussion 
4.1. Phylogeographic Pattern and Processes 


In this study, the lineages A, B, and C of H. chenhsienensis were separated from 
each other over geographic space (Figure 1), and these three lineages were linked by 
54 or 61 mutational steps, demonstrating deep genetic differentiation. Thus, the phylo- 
geographic pattern of H. chenhsienensis conformed to category I defined by Avise [12], 
viz., allopatric lineages with deep differentiation. East Asian monsoon studies in geology 
indicated that the East Asian summer monsoon intensified at about 5.0-4.2 Ma [52-54] 
and 3.7-2.9 Ma [55,56] during the Pliocene (5.333—2.58 Ma), accompanied by noticeably 
increased temperature and monsoon rainfall. These two intensification phases of the East 
Asian summer monsoon in the Pliocene coincided well with the divergence times of three 
lineages for H. chenhsienensis (~4.24 Ma, node 1 in Figure 2; ~3.03 Ma, node 8 in Figure 2). At 
the same time, the global temperature was roughly 2-4 °C warmer in the early mid-Pliocene 
(ca. 5-3 Ma) than at present [57], which resulted in the global sea level being approximately 
10-20 m higher than today [58], suggesting that the basins within distribution regions of H. 
chenhsienensis were isolated from each other during that period. However, H. chenhsienensis 
from each basin did not form a monophyletic group (Figure 1). Hence, two noticeable 
climatic change events in the early mid-Pliocene may be responsible for the divergence of 
lineages A, B, and C and the formation of geographical patterns of H. chenhsienensis. 

The lineage A of H. chenhsienensis consisted of four sub-lineages. Among them, three 
sub-lineages, A-I, A-II, and A-III, were only distributed in the lower reaches of the Yangtze 
River, whereas the sub-lineage A-IV was distributed sporadically in the Ou River, Jiao 
River, and the lower reaches of the Yangtze River (Figure 3a). Paleoclimate studies of East 
Asia showed that the East Asian summer monsoon weakened and the climate became 
cold and dry at about 1.8-1.2 Ma and 0.9-0.3 Ma during the Pleistocene (2.58-0.0117 Ma), 
whereas the East Asian summer monsoon strengthened and the climate turned warm and 
wet at 1.2-0.9 Ma [59]. In addition, all 29 haplotypes of the lineage A were private to 
each sampling site, implicating the low dispersal ability of H. chenhsienensis, inhabiting 
the upper stream. The geographical distribution and the divergence time (1.30, 0.97, and 
0.44 Ma; nodes 2, 4, and 6 in Figure 2, respectively) of the four sub-lineages of lineage 
A in H. chenhsienensis suggested that the low dispersal ability of this species, rather than 
climate change in the Pleistocene, may be the main driving force of sub-lineage splits and 
geographical pattern formation of lineage A. Combining phylogenetic relationships of the 
four sub-lineages of lineage A (Figure 2) and phylogenetic nesting of haplotypes from the 
Ou River and the Jiao River within haplotypes from the lower reaches of the Yangtze River 
(Figure 3a), it was suggested that the lineage A originating from the lower reaches of the 
Yangtze River and H. chenhsienensis from the Ou River and Jiao River in the sub-lineage 
A-IV came from the dispersion of the lower reaches of the Yangtze River population. The 
result of reconstructing palaeodrainages in this study indicated that the Yangtze River, 
Ou River, and Jiao River met the sea independently during the Pleistocene glacial period, 
implying that these three basins were unconnected when the sea level dropped at that 
time. The time of most recent common ancestor of sub-lineage A-IV was 0.28 Ma (node 7 in 
Figure 2), during which the global sea level dropped about 60 m [60], but the climate was 
characterized by warmth and humidity [59,61]. Previous studies have shown that seasonal 
flooding could facilitate the dispersion of fish through a temporary connection between 
rivers that meet the sea independently [62,63]. Consequently, the seasonal floods during 
the late Middle Pleistocene may be the reason for the colonization of H. chenhsienensis from 
the lower reaches of the Yangtze River to the Ou River and Jiao River. 
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The lineage B of H. chenhsienensis was distributed in the Qiantang River and Cao’e 
River (Figure 1). The 18 haplotypes of this lineage were all restricted to each basin, and 
the haplotypes in the two basins were linked by only a single mutational step (Figure 3b), 
which revealed a recent historical connection between the Qiantang River and Cao’e River. 
The TMRCA of lineage B was 0.25 Ma (node 9 in Figure 2), which falls into the late mid- 
Pleistocene glacial period, and the global sea level dropped about 80 m at that time [60]. 
Palaeodrainage reconstruction in this study suggested that the Qiantang River and Cao’e 
River could be interconnected through palaeodrainages as the sea level dropped during the 
glacial period (Figure 1). Therefore, a sea level fall during the late mid-Pleistocene glacial 
period prompted the historical connection between the Qiantang River and Cao’e River, 
and then drove a close relationship between haplotypes of the two basins from lineage B of 
H. chenhsienensis. 

The lineage C of H. chenhsienensis comprises two sub-lineages. The sub-lineage C-I was 
distributed in the Gan River of Poyang Lake drainage (sampling site 20-22), whereas the 
sub-lineage C-II was distributed in the Rao River (sampling site 14), Xin River (sampling 
site 15), and Xiu River (sampling site 16-19) of the Poyang Lake drainage. All 19 haplotypes 
from the two sub-lineages were private haplotypes for each river in the Poyang Lake 
drainage (Table 1), which implies low gene flow among rivers for the lineage C of H. 
chenhsienensis inhabiting the Poyang Lake drainage. The divergence time between sub- 
lineages C-I and C-II was 0.85 Ma (node 10 in Figure 2), when the climate of East Asia 
was cold and dry [59,64]. Consequently, climatic change in the late Early Pleistocene and 
the low dispersal ability of H. chenhsienensis may have driven the lineage divergence and 
geographic pattern of lineage C. 


4.2. Genetic Structure and Historical Demography 


The fixation coefficient of genetic differentiation, Nst, was much higher than Ger, indi- 
cating that the geographical environment was closely related to lineage differentiation in H. 
chenhsienensis. Among basins, high levels of genetic variation detected in pairwise compari- 
son and lack of shared haplotypes reflected limited gene flow in H. chenhsienensis. A high 
genetic differentiation of fish populations from different basins was also observed in the 
Yangtze River and its southern drainage systems, for example, Aphyocypris normalis [14] and 
Hemibagrus guttatus [65]. An increasing number of cases suggested that smaller freshwater 
fish inhabiting the middle and lower water column have poor dispersal ability, which is apt 
to display genetic differentiation, thus resulting in a pronounced genetic structure [66-70]. 
Therefore, the spatial niche for H. chenhsienensis which inhabits the upper reaches of streams 
and the low dispersal ability owing to its small size (<8 cm) may be responsible for the 
strong spatial genetic structure of H. chenhsienensis. 

In this study, the pattern of genetic diversity of H. chenhsienensis in the Yangtze River 
Basin showed high haplotype polymorphism (i) and high nucleotide polymorphism (72), 
which may be ascribed to the simultaneous occurrence of lineages A and C. However, there 
was only a single lineage (A or B) in the Ou River, Qiantang River, or Cao’e River, and 
the pattern of genetic diversity of H. chenhsienensis in each basin was observed as having 
high h values and low 7 values, which suggested that H. chenhsienensis populations in 
these basins may have suffered from genetic bottlenecks, followed by rapid population 
expansions [9]. This speculation was further supported by the borderline significant or sig- 
nificantly negative Fu’s Fs values of lineages A and B, respectively, as well as the unimodal 
peak in mismatch distribution of the lineage B. Moreover, we detected that the population 
expansions of lineages A and B started at ~0.088 Ma or ~0.110 Ma, respectively, but the 
lineage C experienced population contraction from ~0.040 Ma to the present (Figure 4b). 
Paleoclimate studies of East Asia demonstrated that from 0.110 Ma to 0.040 Ma, there was 
an interglacial to glacial transition stage in the Late Pleistocene, when the global sea level 
was 40-50 m lower than the present day, and the climate was cold and dry [71]. Therefore, 
this study concluded that the cross-basin spatial distribution of lineages A and B of H. 
chenhsienensis was the reason for rapid population growth after the lower-temperature 
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period and reduced rainfall during the Late Pleistocene, whereas the spatial distribution of 
the lineage C was confined to the Poyang Lake drainage system, causing a severe popula- 
tion bottleneck after the experience of the dry and cold climate in the Late Pleistocene; thus, 
the population of lineage C had hardly recovered and is in a state of contraction. 


5. Conclusions 


In summary, H. chenhsienensis was composed of three lineages, exhibiting a pattern 
of allopatric lineages with deep divergence. The climate changes during the Pliocene and 
Pleistocene and the poor dispersal ability of H. chenhsienensis might be the major driving 
forces of its phylogeographical pattern. The Yangtze River population of H. chenhsienensis 
denoted high haplotype polymorphism (i) and high nucleotide polymorphism (7z); other 
drainage populations exhibited high h and low 7. During the Late Pleistocene, lineages A 
and B experienced population expansion, whereas lineage C suffered from severe bottleneck 
effects. Based on the spatial subdivision of genetic structure for H. chenhsienensis, it is 
suggested that the four conservation geographic units (Yangtze River, Qiantang River, 
Cao’e River, Jiao River + Ou River) should be classified for management. 
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